NAVAL 

POSTGRADUATE 

SCHOOL 

MONTEREY,  CALIFORNIA 


THESIS 


ANGULAR  RATE  ESTIMATION  BY 
MULTIPLICATIVE  KALMAN  FILTERING 
TECHNIQUES 

by 

Vincent  C.  Watson 
December  2003 

Thesis  Advisor:  Roberto  Cristi 

Co-Advisor:  Brij  Agrawal 


Approved  for  Public  Release;  Distribution  is  Unlimited 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  No.  0704-0188 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including 
the  time  for  reviewing  instruction,  searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and 
completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any 
other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington 
headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite 
1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project 
(0704-0188)  Washington  DC  20503,  _ _ 

I.  AGENCY  USE  ONLY  (Leave  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

December  2003  Master’s  Thesis 

4.  TITLE  AND  SUBTITLE:  Mitigation  of  Spacecraft  Attitude  Estimation  5.  FUNDING  NUMBERS 
Error  Via  Kalman  Filtering 

6.  AUTHOR(S)  Vincent  C.  Watson  _ 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES)  8.  PERFORMING 

Naval  Postgraduate  School  ORGANIZATION  REPORT 

Monterey,  CA  93943-5000  NUMBER 

9.  SPONSORING  /MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES)  10.  SPONSORING/MONITORING 
N/A  AGENCY  REPORT  NUMBER 

II.  SUPPLEMENTARY  NOTES  The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official 

policy  or  position  of  the  Department  of  Defense  or  the  U.S.  Government. _ 

12a.  DISTRIBUTION  /  AVAILABILITY  STATEMENT  12b.  DISTRIBUTION  CODE 

Approved  for  Public  Release;  Distribution  is  Unlimited 

13.  ABSTRACT  (maximum  200  words) 

Spacecraft  attitude  estimation  and  pointing  accuracy  have  always  been  limited  by  imperfect  sensors.  The 
rate  gyroscope  is  one  of  the  most  critical  instruments  used  in  spacecraft  attitude  estimation  and  unfortunately 
historical  trends  show  this  instrument  degrades  significantly  with  time.  Degraded  rate  gyroscopes  have  impacted 
the  missions  for  several  NASA  and  ESA  spacecraft,  including  the  Hubble  Telescope.  A  possible  solution  to  this 
problem  is  using  a  mathematically  modeled  dynamic  gyroscope  in  lieu  of  a  real  one.  In  this  thesis,  data  from  such 
a  gyro  is  presented  and  integrated  into  a  spacecraft  attitude  estimation  algorithm. 

The  impediment  to  spacecraft  attitude  estimation  presented  by  imperfect  sensors  has  been  overcome  by 
developing  more  accurate  sensors  and  using  Kalman  filters  to  reduce  the  effect  of  noisy  measurements.  Kalman 
filters  for  spacecraft  attitude  estimation  have  historically  been  based  on  an  Euler  angle  or  quaternion  formulation. 
Though  Euler  angles  and  quaternions  are  arguably  the  easiest  methods  with  which  to  describe  the  attitude  of  a 
spacecraft,  other  methods  of  describing  attitudes  do  exist  -  including  the  Gibbs  and  Rodriguez  parameters.  A 
Kalman  filter  based  upon  the  Gibbs  parameter  is  presented  and  analyzed  in  this  thesis. 


14.  SUBJECT  TERMS  Kalman  Filter,  Gibbs  Parameter,  Dynamic  Gyroscope,  Attitude  Estimation  15.  NUMBER  OF 

PAGES 
_ 70 

16.  PRICE  CODE 


17.  SECURITY 

18.  SECURITY 

19.  SECURITY 

20.  LIMITATION 

CLASSIFICATION  OF 

CLASSIFICATION  OF  THIS 

CLASSIFICATION  OF 

OF  ABSTRACT 

REPORT 

PAGE 

ABSTRACT 

Unclassified 

Unclassified 

Unclassified 

UL 

NSN  7540-0 1  -280-5500  Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Std.  239-18 


1 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


11 


Approved  for  Public  Release;  Distribution  is  Unlimited 


ANGULAR  RATE  ESTIMATION  BY  MULTIPLICATIVE  KALMAN 
FILTERING  TECHNIQUES 

Vincent  C.  Watson 
Lieutenant,  United  States  Navy 
B.S.,  Georgia  Institute  of  Technology,  1994 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF  SCIENCE  IN  ASTRONAUTICAL  ENGINEERING 


from  the 


NAVAL  POSTGRADUATE  SCHOOL 
December  2003 


Author:  Vincent  C.  Watson 


Approved  by:  Roberto  Cristi 

Thesis  Advisor 


Brij  N.  Agrawal 
Co-Advisor 


Anthony  J.  Healey 

Chairman,  Department  of  Mechanical  and  Astronautical 
Engineering 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


IV 


ABSTRACT 


Spacecraft  attitude  estimation  and  pointing  accuracy  have  always  been  limited  by 
imperfect  sensors.  The  rate  gyroscope  is  one  of  the  most  critical  instruments  used  in 
spacecraft  attitude  estimation  and  unfortunately  historical  trends  show  this  instrument 
degrades  significantly  with  time.  Degraded  rate  gyroscopes  have  impacted  the  missions 
for  several  NASA  and  ESA  spacecraft,  including  the  Hubble  Telescope.  A  possible 
solution  to  this  problem  is  using  a  mathematically  modeled  dynamic  gyroscope  in  lieu  of 
a  real  one.  In  this  thesis,  data  from  such  a  gyro  is  presented  and  integrated  into  a 
spacecraft  attitude  estimation  algorithm. 

The  impediment  to  spacecraft  attitude  estimation  presented  by  imperfect  sensors 
has  been  overcome  by  developing  more  accurate  sensors  and  using  Kalman  filters  to 
reduce  the  effect  of  noisy  measurements.  Kalman  filters  for  spacecraft  attitude 
estimation  have  historically  been  based  on  an  Euler  angle  or  quaternion  formulation. 
Though  Euler  angles  and  quaternions  are  arguably  the  most  common  methods  with  which 
to  describe  the  attitude  of  a  spacecraft,  other  methods  of  describing  attitudes  do  exist  - 
including  the  Gibbs  and  Rodriguez  parameters.  A  Kalman  filter  based  upon  the  Gibbs 
parameter  is  presented  and  analyzed  in  this  thesis. 
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I.  INTRODUCTION 


A.  SPACECRAFT  ATTITUDE  ESTIMATION 

The  pointing  accuracy  requirements  of  modem  Department  of  Defense  satellites 
have  increased  steadily  over  the  past  few  decades.  As  the  pointing  accuracy  requirements 
for  a  satellite  increase,  so  does  the  level  of  accuracy  in  the  estimation  of  its  attitude. 
Spacecraft  attitude  estimation  is  an  extremely  complex  and  non-linear  process.  Inputs 
from  several  different  sensors  -  such  as  star  trackers  and  rate  gyroscopes,  among  others  - 
are  required  to  determine  the  attitude  of  a  spacecraft.  Because  no  sensor  is  error-free  - 
even  at  the  time  of  manufacture  -  and  all  man-made  equipment  degrades  with  age,  the 
problem  of  accurate  attitude  estimation  throughout  the  mission  life  of  a  satellite  is 
extremely  important.  A  method  for  overcoming  these  impediments  via  a  Kalman 
filtering  process  is  presented  and  analyzed  here. 

Kalman  filtering  has  been  used  in  spacecraft  attitude  estimation  for  quite  some 
time.  The  earliest  published  application  was  in  1970  by  Farrell  and  several  others  have 
followed  since.  Lefferts,  Shuster,  and  Markley  published  a  thorough  review  of  the  topic 
in  1982  and  since  then  Markley,  Crassidis,  and  several  others  have  kept  Kalman  filtering 
an  active  topic  of  research  in  the  space  industry  [Markley].  Several  different  attitude 
representations  have  been  used  in  Kalman  filtering  with  varying  degrees  of  success. 
Palenno  successfully  implemented  a  Kalman  filter  using  an  Euler  angle  representation  of 
the  attitude  for  a  simulated  bifocal  relay  mirror  spacecraft  [Palermo]  and  his  dynamic 
model  is  used  as  a  starting  point  for  this  work. 

B.  SATELLITE  ENVIRONMENT 

The  environment  in  which  a  satellite  operates  complicates  the  problem  of 
spacecraft  attitude  estimation  and  control.  A  satellite  in  orbit  is  subject  to  non-constant 
external  torques  at  all  times  -  some  secular  (varying  linearly  with  time),  some  periodic, 
and  some  random.  The  forces  that  induce  these  torques  are  gravity  gradient,  magnetic, 
solar  pressure,  and  atmospheric  drag,  though  the  effects  of  atmospheric  drag  on  satellites 
outside  of  LEO  is  considered  negligible.  Though  these  torques  can  be  predicted  with 
high  accuracy  (with  the  exception  of  atmospheric  drag),  any  errors  in  prediction  couple 
with  errors  in  sensor  accuracy  to  increase  the  error  in  attitude  estimation.  Because  large 
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errors  in  attitude  can  be  attributed  to  the  effects  of  these  external  torques,  a  brief 
overview  of  gravity  gradient,  magnetic,  and  solar  torque  and  the  errors  that  can  result 
from  each  is  presented  here.  For  more  in  depth  information  the  reader  is  referred  to  texts 
by  Sidi,  Bong  Wie,  or  Hughes  [Sidi,  Bong  Wie,  Hughes].  In  the  rest  of  this  section  we 
survey  a  number  of  disturbances  acting  on  the  spacecraft. 

1.  Gravity  Gradient  Torque 

A  spacecraft  is  not  a  point  mass  and  may  not  be  treated  as  such.  It  is  a  rigid  body 
(in  most  cases)  with  a  mass  distribution  about  a  center  of  mass.  Vice  treating  it  as  a  point 
mass  an  inertia  dyadic  is  used.  Gravity  gradient  torques  are  imparted  to  a  spacecraft 
because  gravity  acts  on  each  element  of  the  spacecraft.  Gravity  acting  on  a  mass  m 
located  at  a  distance  r  from  the  spacecraft  center  of  mass  will  induce  a  torque  about  the 
center  of  mass.  The  effects  of  gravity  will  act  on  each  portion  of  the  satellite  in 
accordance  with  Newton’s  Laws  of  Gravitation.  Ignoring  the  effects  of  the  moon  and 
other  third  bodies  due  to  their  small  effects,  these  torques  may  be  written  as: 

A  =  (i-i) 

where  //  is  the  gravitational  parameter  for  the  Earth,  R(l  is  the  distance  from  the  satellite 
to  the  center  of  the  Earth,  and  I B  is  the  inertia  dyadic  for  the  satellite. 

Some  satellites  with  low  pointing  requirements  actually  use  gravity  gradient 
stabilization.  A  quick  examination  of  equation  1-1  reveals  that  knowledge  of  the  position 
of  the  satellite  and  the  inertia  dyadic  are  critical  for  predicting  the  gravitational  torque. 
System  identification  algorithms  are  being  actively  researched  and  developed  to  facilitate 
more  precise  knowledge  of  spacecraft  inertia  dyadics. 

In  addition  to  its  own  mass  distribution,  the  fact  that  the  gravitational  field 
produced  by  the  Earth  is  an  aspherical  potential  further  complicates  the  problem.  Like  a 
satellite,  the  Earth  does  not  have  a  uniform  mass  distribution  -  hence  its  aspherical 
gravity  field.  Zonal,  tesseral,  and  sectoral  harmonics  within  the  field  are  both  non-linear 
and  extremely  complex  to  model  mathematically.  Inaccuracies  in  modeling  the 
gravitational  field  of  the  Earth  can  lead  to  attitude  pointing  errors  as  well. 
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2.  Magnetic  Torque 

The  Earth  has  a  rotating  magnetic  field  B  through  which  any  orbiting  satellite 


must  travel.  The  satellite  itself  has  an  intrinsic  magnetic  moment  in  = 


mx 

m 

m. 


.  Interaction 


between  the  magnetic  moment  of  the  spacecraft  and  the  magnetic  field  of  the  Earth 
generates  a  torque  on  the  spacecraft  calculated  by 


Nb  =  in  xB 


(1-2) 


The  magnitude  of  B  will  decrease  as  the  satellite  altitude  increases.  For  dealing 
with  torques  induced  by  the  interaction  of  the  spacecraft  and  Earth  magnetic  fields, 
precise  knowledge  of  the  spacecraft  magnetic  moment  m  is  critical.  As  with  the  moment 
of  inertia,  system  identification  algorithms  can  be  used  to  mitigate  errors  due  to 
inaccurate  estimates  of  m  . 

3.  Solar  Torque 

Maxwell’s  equations  imply  that  electromagnetic  waves  have  momentum,  which 
may  be  transferred  to  objects  with  which  it  comes  in  contact.  Since  light  is  an 
electromagnetic  wave,  it  exerts  pressure.  Though  this  pressure  is  miniscule  in  an  Earth 
environment  it  is  not  miniscule  for  a  satellite  in  orbit  about  the  Earth. 

For  electromagnetic  radiation,  basic  physics  shows  us  that  Electromagnetic  Force 
=  Work  +  Energy  Density  which,  when  expressed  in  a  more  mathematical  manner 
becomes 

[  (TUda)da  =  [  FdV  +  —  f  —^dV  (1-3) 

is  iv  dtJl'c~ 

where  T  is  the  Maxwell  stress  tensor,  F  is  the  force  density,  and  S  is  the  Poynting 

vector  (S  =  —ExB ). 

An 

The  amount  of  pressure  exerted  on  an  object  by  an  electromagnetic  wave  is  highly 
dependent  upon  the  type  of  surface  being  illuminated.  In  some  complex  models,  the  type 
of  reflection  -  specular  or  diffuse  -  as  well  as  the  reflectivity  of  the  surface  plays  a  part  in 
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the  calculation.  In  a  simple  model  the  following  equation  is  often  used  to  determine  the 
maximum  solar  torque  imparted  to  a  spacecraft 


N, 


SP 


F 

=  —As(\  +  q) cos(/)(c  - c  ) 
c 


(1-4) 


where  As  is  the  spacecraft  area,  q  is  the  reflectance  factor,  i  is  the  angle  of 
incidence  of  the  incoming  light,  c  -c  is  the  distance  between  the  spacecraft  center  of 
gravity  and  the  center  of  solar  pressure,  and  Fs  is  the  solar  power  density  (which  varies 
with  time).  Solar  cycles,  changes  in  csp  -c  due  to  fuel  expenditures,  and  the  changing 

area  of  a  spacecraft  tracking  a  point  on  the  Earth  all  contribute  to  the  difficulty  of 
estimating  and  compensating  for  disturbances  due  to  solar  pressure  [SMAD]. 

4.  Atmospheric  Drag 

As  previously  mentioned  the  effects  of  atmospheric  drag  on  satellites  outside  of 
LEO  are  considered  negligible.  For  LEO  satellites,  however,  atmospheric  drag  is  the 
most  difficult  external  torque  to  predict.  Because  the  dynamics  of  the  outer  reaches  of 
the  atmosphere  are  not  fully  understood  it  cannot  be  modeled  accurately.  The  effects  of 
drag  on  LEO  spacecraft  are  directly  proportional  to  the  area  of  the  spacecraft. 
Atmospheric  drag,  though  pertinent  to  LEO  applications,  was  not  included  in  any  of  the 
models  used  in  developing  this  thesis  but  is  mentioned  here  for  completeness. 

C.  SENSORS 

Accurate  measurement  of  data  from  external  sources  is  required  for  a  spacecraft 
attitude  control  system  to  estimate  its  attitude.  Though  several  sensors  exist  that  perform 
this  function  -  including  sun  sensors,  horizon  sensors,  and  earth  sensors  among  others  - 
the  star  tracker  is  the  most  accurate  and  pertinent  to  the  focus  of  this  work. 

1.  Rate  Gyroscopes 

A  cursory  perusal  of  either  the  Euler  Equations  of  Motion  or  the  kinematic 
equations  for  the  quaternion  or  Gibbs  parameters  reveals  their  dependence  upon  angular 
rate  data.  Any  error  in  angular  rate  measurement  or  calculation  will  result  in  an  error  in 
spacecraft  attitude  estimation.  Angular  rate  information  is  critical  to  the  accurate 
estimation  of  spacecraft  attitude  -  regardless  of  the  estimation  method.  Rate  gyroscopes 
provide  this  information  to  the  spacecraft  attitude  control  system. 
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There  are  several  different  kinds  of  rate  gyroscopes  available  for  use  in  the  space 
environment  today.  Mechanical  rate  gyroscopes  and  rate  integrating  gyroscopes  have 
been  in  use  on  orbit  for  quite  several  years.  Both  of  these  types  of  gyroscopes  are 
dependent  upon  mechanical  moving  parts  which  degrade  over  time.  Laser  gyroscopes, 
quartz  rate  sensors,  and  hemispherical  resonator  gyroscopes  are  now  available  today  with 
much  higher  accuracies  and  no  reliance  on  moving  parts  [Sidi]. 

2.  Dynamic  Gyroscopes 

Rate  gyroscopes  are  man-made  devices  and  therefore  have  finite  lifetimes  - 
especially  in  the  harsh  environs  of  space.  Because  all  of  the  external  secular  and  periodic 
torques  on  a  spacecraft  may  be  modeled  with  some  modicum  of  accuracy  and  the  torques 
applied  to  the  spacecraft  via  the  momentum  exchange  devices  of  its  attitude  control 
system  are  known  with  a  fair  amount  of  precision,  it  is  possible  to  detennine  the  angular 
rate  of  the  spacecraft  via  mathematical  modeling  and  the  use  of  external  measurements. 
Such  an  algorithm  is  called  a  dynamic  gyroscope.  The  use  of  a  dynamic  gyroscope  upon 
failure  of  a  simulated  mechanical  gyroscope  is  shown  and  analyzed  in  this  work. 

3.  Star  Trackers 

Navigators  from  ancient  times  used  the  stars  as  navigational  aids.  Satellites 
navigating  in  space  do  the  same  via  a  device  called  a  star  tracker.  Because  stars  may  be 
considered  inertially  fixed  bodies  for  all  intents  and  purposes  and  because  they  are 
extremely  small  as  seen  from  our  solar  system,  they  are  ideal  objects  to  use  as  an  attitude 
reference.  Star  trackers  provide  the  most  accurate  attitude  data  to  the  spacecraft  attitude 
control  system  by  several  orders  of  magnitude  over  any  other  type  of  sensor.  While  older 
star  trackers  were  capable  of  tracking  only  one  star  at  a  time,  the  new  generation  of  star 
trackers  can  feed  attitude  quaternions  to  its  host  satellite  and  track  multiple  stars 
simultaneously.  A  new  type  of  star  sensor  is  currently  being  developed  by  Dr.  Junkins  et 
al  which  will  use  two  simultaneous  images  of  star  fields  from  different  parts  of  the  sky  to 
detennine  spacecraft  attitude  [Junkins]. 

D.  KALMAN  FILTERS 

As  previously  mentioned,  a  satellite  attitude  control  system  receives  input  from 
imperfect  external  sensors  to  estimate  the  actual  attitude  of  the  satellite.  The  Kalman 
filter  is  an  estimation  algorithm  frequently  employed  to  do  this.  In  essence,  the  Kalman 
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filter  is  a  set  of  mathematical  equations  that  provides  a  recursive  solution  of  the  least- 
squares  estimation  method.  It  supports  computational  estimations  of  past,  present,  and 
future  states  and  can  do  so  even  when  the  nature  of  the  system  under  consideration  is  not 
precisely  known. 

Given  the  extremely  non-linear  and  varying  nature  of  the  space  environment,  this 
makes  it  an  ideal  tool  with  which  to  perfonn  spacecraft  attitude  estimation.  The 
effectiveness  of  the  filtering  algorithm  is  highly  dependent  upon  how  the  state  vector  is 
defined,  the  type  of  filter  employed,  and  several  other  factors.  In  this  research  we 
developed  a  Kalman  Filter  based  on  specific  parameterization  of  the  spacecraft  attitude 
and  attitude  error. 

E.  STAR  GAPS  AND  RATE  GYROSCOPE  UPSETS 

Satellites  in  orbit  currently  experience  periods  in  which  their  star  trackers  are 
unable  to  sense  stars  -  called  star  gaps.  Rate  gyroscopes  in  some  older  satellites  are 
sometimes  sending  data  over  1000  times  the  actual  reading  -  such  an  event  is  called  a 
rate  gyroscope  upset.  Both  of  these  events  wreak  havoc  upon  attitude  control  algorithms 
-  the  Kalman  filter  in  particular.  Both  of  these  phenomena  are  simulated  and  their  effects 
mitigated  via  adaptive  covariance  and  dynamic  gyroscope  integration. 

F.  THESIS  ORGANIZATION 

In  this  thesis  we  introduce  the  dynamics  of  spacecraft  attitude  in  Chapter  Two. 
Chapter  Three  begins  with  a  cursory  overview  of  the  Kalman  filtering  process  followed 
by  a  brief  section  on  attitude  error  representation.  It  closes  with  a  derivation  of  the  Gibbs 
parameter  based  Kalman  filter.  Chapter  Four  covers  the  problems  caused  by  star  gaps 
and  methods  of  managing  such  problems.  Chapter  Five  introduces  the  dynamic 
gyroscope  and  methods  of  dealing  with  rate  gyroscope  upsets. 
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II.  SPACECRAFT  ATTITUDE 


A.  ATTITUDE  PARAMATRIZATON 

In  this  section  we  review  spacecraft  attitude  estimation  methods.  There  are 
several  methods  for  doing  this,  but  only  a  short  description  of  most  of  them  will  be 
included  here.  Emphasis  will  be  placed  on  the  quaternion  because  understanding  the 
quaternion  is  critical  to  comprehending  the  development  of  the  Kalman  filters  used  in 
spacecraft  attitude  estimation. 

Regardless  of  the  method,  describing  the  attitude  of  a  spacecraft  basically  entails 
describing  the  orientation  of  one  coordinate  system  with  respect  to  another.  For  purposes 
of  satellites,  a  coordinate  system  based  on  the  principal  axes  of  the  spacecraft  -  called  the 
body  coordinate  system  -  and  Earth  Centered  Inertial  (ECI)  systems  are  used. 

1.  Direction  Cosine  Matrix 

The  direction  cosine  matrix  is  the  simplest  manner  in  which  to  describe  the 
attitude  of  a  spacecraft.  Given  two  coordinate  systems  each  consisting  of  three 
orthogonal  unit  vectors  there  exists  a  three  by  three  matrix  C  that  relates  the  two 
coordinate  systems. 

When  dealing  with  different  coordinate  systems  it  becomes  necessary  to  develop 
a  notation  for  annotating  in  which  coordinate  system  a  vector  is  expressed.  In  this  paper 
a  superscript  to  the  left  of  the  vector  will  indicate  the  coordinate  system  in  which  a  vector 
is  expressed.  For  example,  'v  is  a  vector  in  the  inertial  reference  frame  while  hv  is  a 
vector  in  the  body  reference  frame. 


A  vector  x  = 


expressed  in  body  coordinates  may  be  expressed  in  inertial 


coordinates  by  multiplying  it  by  a  direction  cosine  matrix  (DCM)  as  shown  below 

iCbbx  =  ix  (2-1) 


Note  the  two  superscripts  on  the  DCM.  For  a  direction  cosine  matrix,  the 
superscript  on  the  left  indicates  the  coordinate  system  the  DCM  transfonns  a  vector  to 
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while  the  right  superscript  indicates  the  coordinate  system  from  which  the  vector  is  being 
transformed.  Conversely,  a  vector  in  the  inertial  frame  may  be  expressed  in  body  frame 
coordinates  by 

bC‘  ‘x  =  hx  (2-2) 

As  one  might  expect  'Ch  and  hC‘  are  related  as  the  transpose  of  each  other. 

2.  Euler  Angles 

Euler  angles  are  a  series  of  rotations  about  the  axes  of  one  coordinate  system  to 
align  it  to  another  coordinate  system.  These  rotations  are  sometimes  referred  to  as  roll  (p 
(rotation  about  the  x  axis),  pitch  0  (rotation  about  the  y  axis),  and  yaw  i//  (rotation  about 
the  z  axis). 

Since  in  the  Euler  angle  formulation  the  order  matters,  there  are  twelve  possible 
combinations  of  Euler  angles  one  may  use  to  go  from  one  coordinate  system  to  another. 
Multiplication  of  the  individual  rotation  matrices  for  each  Euler  angle  will  result  in  a 
direction  cosine  matrix.  Given  the  two  coordinate  systems  a  and  b  the  direction  cosine 
matrix  for  a  roll-pitch-yaw  Euler  sequence  would  be 

1  0  0  cos  0  0  -sin  0  cos y  sin  y/  0 

aCh  =  0  cos  cp  sin  <p  0  1  0  -sin^  cos  y/  0  (2-3) 

0  -sin  cp  cos  cp  sin6*  0  cos6*  0  0  1 

3.  Quaternions 

A  mathematical  structure  called  the  quaternion  is  a  convenient  method  with  which 
to  describe  the  orientation  of  a  coordinate  system. 

a.  Quaternion  Definition 

Quaternions  are  a  fonn  of  hyper-complex  numbers  invented  by  William 
Hamilton  in  the  nineteenth  century  and  are  today  used  extensively  both  in  robotics  and 
spacecraft  attitude  control.  Quaternions  are  represented  as  a  four  element  set  consisting 
of  three  vector  components  and  one  scalar  component 
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There  is  no  set  convention  for  which  element  of  the  quaternion  is  a  scalar, 
but  in  this  paper  q4  will  always  be  the  scalar  element  of  the  quaternion. 
b.  Attitude  Quaternions 

When  used  for  describing  an  attitude,  quaternions  are  constrained  to  the 
surface  of  a  four  dimensional  hyper-sphere  defined  by 

qf  +  q;  +  q]  +  ql=\  (2-5) 

To  determine  the  physical  meaning  of  quaternions  constrained  on  this 
hyper-sphere,  consider  two  coordinate  systems  as  shown  below,  S  and  S’.  The  vector 
portion  of  the  quaternion  represents  an  axis  about  which  coordinate  system  S  must  be 
rotated  to  align  it  with  the  5”  coordinate  system  -  called  the  eigenaxis. 

s 


Figure  1.  Coordinate  Systems  S  and  S’ 


The  scalar  portion  of  the  quaternion  is  a  measure  of  the  magnitude  of  the 


rotation#,. .  Defining  the  eigenaxis  as  a  unit  vector?  = 


be  defined  as  follows: 


ql  =  sin 
q2  =  sin 
q3  =  sin 
q4  =  cos 


(oA 
l  2  J 
(oA 

l  2  J' 

1 2  j 


,  the  quaternion  elements  may 


(2-6) 


Since  quaternions  are  vectors  of  unit  magnitude,  they  may  not  be  added 
together  through  the  standard  definition  of  addition.  Quaternion  algebra  is  a  rich  and 
interesting  topic,  but  one  that  will  not  be  dealt  with  here  in  its  entirety.  Two  uses  of 
quaternion  algebra  will  be  addressed  because  of  their  usefulness  in  attitude  estimation  - 
the  quaternion  conjugate  and  quaternion  multiplication. 

c.  The  Quaternion  Conjugate 

The  quaternion  conjugate  of  quaternion  q  is  given  by 

q  =  -<h-q2  ~<h  +  q*  (2-2) 


From  a  physical  standpoint,  the  conjugate  of  a  quaternion  q  represents  a 
rotation  of  the  same  magnitude  about  a  vector  in  the  opposite  direction.  It  is  also  worthy 
to  note  that  the  inverse  of  a  quaternion  is  identical  to  its  complex  conjugate.  The 
quaternion  conjugate  is  useful  when  using  the  quaternion  as  a  rotation  operator,  which 
will  be  discussed  later. 

d.  Quaternion  Multiplication 

As  noted  previously,  quaternions  have  their  own  unique  algebra.  The 
symbol  0  is  used  to  denote  quaternion  multiplication.  It  can  be  shown  that  the 
quaternion  product  is  defined  by 
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p®q  = 


d4 

<h 

~q2 

<h 

Pi 

-q3 

di 

q2 

Pi 

<h 

-9, 

d4 

<h 

Pi 

-q2 

-<h 

q4_ 

_P4_ 

(2-8) 


which  may  also  be  written  as 


p<$q  = 


p4q+q4p~pxq 

p4q4-p~q 


(2-9) 


The  quaternion  product  is  a  useful  tool  for  time  propagation  of  the 
quaternion  -  as  will  be  shown  later.  Recall  that  a  quaternion  represents  a  rotation  from  a 
reference  frame  to  a  given  attitude.  It  naturally  follows  that  the  quaternion  product 
defined  in  equation  2-9  would  also  represent  a  rotation  -  which  it  does.  Defining  A{q)  as 
a  DCM  representing  the  same  rotation  as  quaternion  q,  it  may  be  shown  that 

A(p®q)  =  A(p)A(q)  (2-10) 

e.  Quaternion  Propagation 

Given  an  angular  velocity  vector  d>  in  body  coordinates  and  an  initial 
orientation  expressed  by  the  quaternion  q0  the  orientation  at  any  time  t  may  be  expressed 
as  the  solution  of  the  following  differential  equation 


dq 

dt 


=  q  = 


(2-11) 


with  the  initial  condition  q(0)  =  q0 . 

f  Advantages  of  Quaternions 

Quaternion  representations  of  spacecraft  attitude  hold  several  advantages 
over  the  Euler  angle  and  direction  cosine  matrix  representations.  Perhaps  the  most 
obvious  advantage  comes  from  the  size  -  a  quaternion  is  a  four  element  structure, 
whereas  a  direction  cosine  matrix  has  nine  elements.  A  quaternion  describes  the  same 
attitude  with  half  of  the  number  of  elements,  saving  both  memory  and  processing  power. 
The  simplicity  of  the  time  derivative  of  the  quaternion  makes  it  an  ideal  method  with 
which  to  do  attitude  propagation.  Also  of  great  importance  is  the  fact  that  it  has  no 
singularities. 
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4.  Gibbs  Vector 

Gibbs  vectors  are  closely  related  to  quaternions,  though  they  exist  in  □  3  vice  4 . 
The  Gibbs  vector  is  defined  as 


_  1 

g  =  — 

^4 


^2 


=  e  tan 


(<C\ 

UJ 


(2-12) 


The  Gibbs  vector  is  obviously  singular  whenever  q4  goes  to  zero,  which  occurs 
for  rotations  of  180  degrees.  Consequently,  the  Gibbs  parameter  is  an  extremely 
effective  tool  for  describing  rotations  in  the  interval  (-180°, +  180°) .  The  relationship 
between  the  Gibbs  vector  and  the  quaternion  is  shown  in  the  gnomonic  projection  in 
figure  2: 


Figure  2.  Gibbs  vector  as  a  Gnomonic  Projection 

B.  ATTITUDE  SENSORS 

Measurement  of  external  data  is  required  in  order  to  produce  an  attitude  estimate. 
There  are  a  multitude  of  sensors  available  to  measure  attitude  data  for  a  spacecraft, 
including  horizon  sensors,  sun  sensors  (coarse  and  fine),  and  star  trackers.  None  of  these 
provides  a  perfect  measurement,  hence  the  need  for  attitude  error  estimation.  A  brief 
overview  of  rate  gyroscopes  and  star  trackers  along  with  their  respective  sources  of  error 
is  presented  here,  as  both  are  simulated  in  models  used  in  this  thesis. 


12 


1.  Rate  Gyroscopes 

Rate  gyroscopes  provide  angular  rate  data  to  the  spacecraft  attitude  control 
system.  There  are  several  different  types  of  rate  gyroscopes  available  today,  including 
mechanical,  fiber  optic,  laser,  and  hemispherical  resonator  gyroscopes.  Regardless  of  the 
type  of  gyroscope,  all  perfonn  the  same  function.  Each  sends  the  spacecraft  attitude 
control  system  the  angular  rates  of  the  satellite  in  body  coordinates.  It  is  beyond  the 
scope  of  this  work  to  describe  each  in  detail,  but  one  thing  that  all  rate  gyroscopes  have  in 
common  is  that  they  are  inherently  noisy.  The  errors  sent  to  the  attitude  control  system 
by  noisy  rate  gyroscopes  may  be  categorized  into  four  distinct  types  of  error:  bias  error, 
random  walk,  scale  factor  error,  and  orthagonality  errors. 

a.  Bias 

Bias  error  is  sometimes  referred  to  as  a  constant  drift.  Bias  errors  for 
current  generation  gyroscopes  range  from  0.00  T  /  hour  to  T  /  hour .  The  measured 
angular  rate  for  a  rate  gyroscope  may  be  expressed  as 

(°meaS(t)  =  ®(t)  +  b(t )  (2-13) 

With  b(t)  representing  the  bias  at  time  t.  As  will  be  shown  later  in  this 
work,  it  is  possible  to  estimate  gyroscope  bias  via  a  Kalman  filtering  process. 

b.  Rate  Walk 

Rate  walk  is  sometimes  described  as  bias  drift.  Random  walk  propagates 
from  the  white  noise  of  the  sensor  and  is  responsible  for  non-detenninistic  behavior.  Its 

units  are  normally0/ yjhour  . 

c.  Scale  Factor 

The  scale  factor  error  is  the  linear  deviation  of  the  measured  rate  from  the 
true  rate  -  normally  given  as  a  percentage  or  in  parts  per  million.  Asymmetry  and  non¬ 
linearity  have  been  observed  in  scale  factors  for  rate  gyroscopes.  Scale  factors  are  caused 
by  imperfections  in  manufacturing  and  the  degradation  of  the  gyroscope  with  the  passage 
of  time  [Hewitson  et  al].  Though  a  source  of  error  in  any  gyroscope,  scale  factors  are  not 
used  in  any  of  the  models  in  this  work,  but  are  presented  here  for  completeness. 
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d.  Orthogonality  Errors 

Rate  gyroscopes  are  mounted  in  clusters  relative  to  one  another  to  provide 
data  to  a  spacecraft  attitude  control  system.  Any  misalignment  of  the  sensors  with 
respect  to  one  another  will  result  in  error  -  known  as  orthogonality  error.  Launch 
vibrations  and  thermal  deformation  over  time  could  cause  such  errors  to  occur.  As  with 
scale  factor  errors,  no  orthogonality  errors  are  modeled  in  this  work  -  they  are  presented 
for  theoretical  completeness. 

2.  Star  Trackers 

As  their  name  implies,  star  trackers  generate  attitude  data  based  upon  the  relative 
position  of  stars  as  seen  from  the  spacecraft  itself.  Each  star  sensor  has  a  star  catalog 
with  which  to  compare  the  image  it  sees.  Due  to  the  large  number  of  recorded  stars,  the 
exact  composition  of  the  catalog  will  be  determined  by  the  orbit  in  which  the  host 
satellite  will  be  placed  -  there  is  no  such  thing  as  a  standard  star  catalog. 

a.  Star  Characteristics 

In  order  to  use  stars  for  navigational  purposes  it  is  necessary  to 
differentiate  between  stars.  There  are  two  characteristics  of  stars  which  make  this 
possible  -  magnitude  and  spectra. 

Magnitude  refers  to  the  brightness  of  a  star.  There  are  two  types  of 
magnitude  -  absolute  and  apparent.  The  apparent  magnitude  is  how  bright  an  object 
seems  when  viewed  from  the  Earth.  In  contrast,  the  absolute  magnitude  is  the  apparent 
magnitude  of  an  object  placed  at  ten  parsecs  away.  Star  trackers  use  the  apparent 
magnitude  for  their  calculations. 

The  spectrum  of  a  star  refers  to  the  type  of  radiation  it  is  emitting.  Stellar 
spectra  are  divided  into  seven  categories.  There  are  O,  B,  A,  F,  G,  K,  and  M  class  stars, 
with  O  being  the  hottest  and  M  being  the  coolest  star. 

b.  Star  Sensing 

Current  generation  star  trackers  use  Charge  Coupled  Devices  (CCDs)  to 
detect  stars,  though  future  generation  star  trackers  will  use  CMOS  technology  [Junkins]. 
A  stray  light  shield  prevents  light  from  outside  the  star  tracker  bore  sight  from  reaching 
the  CCD.  Once  the  CCD  has  registered  a  star  field  image,  the  star  sensor  processor  will 
match  that  image  with  its  library  and  send  an  attitude  quaternion  to  the  spacecraft  attitude 
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control  system.  This  is  an  extremely  simplified  explanation  of  what  actually  occurs,  but 
the  actual  mechanics  of  this  process  are  outside  the  scope  of  this  work. 

c.  Star  Tracking 

Once  one  or  more  stars  have  been  acquired,  the  spacecraft  will  move 
around  the  Earth  in  its  orbit  and  possibly  change  its  attitude  intentionally.  The  software 
must  track  the  image  of  the  star  through  this  motion  and  be  ready  to  acquire  a  new  star 
when  the  image  exits  the  Field  of  View  (FOV)  of  the  star  tracker  [Sidi] . 

d.  Star  Gaps 

There  are  times  when  a  star  exits  the  field  of  view  (FOV)  of  a  star  tracker 
and  no  new  star  within  the  new  FOV  has  been  acquired.  During  such  time  periods  no 
data  is  being  sent  to  the  satellite  attitude  control  system  -  such  periods  are  called  star 
gaps.  Star  gaps  may  also  be  caused  by  the  satellite  pointing  in  such  a  way  that  the  moon 
or  sun  prevents  the  star  tracker  from  sensing  any  stars  due  to  their  high  intensity.  Star 
gaps  are  detrimental  to  attitude  detennination  when  they  extend  for  more  than  a  few 
seconds  and  can  wreak  havoc  on  attitude  estimation  algorithms,  as  will  be  shown  in  a 
later  section. 

e.  Star  Tracker  Error 

Although  easily  the  most  accurate  sensor  available,  star  trackers  are  not 
error  free.  Fike  rate  gyroscopes,  star  trackers  also  have  a  scale  factor,  which  will  vary 
with  age.  The  CCD  assemblies  are  temperature  sensitive  structures  -  a  primary  source  of 
error.  Misalignments  of  star  trackers  with  respect  to  the  spacecraft  body  as  a  result  of 
launch  vibration  or  thennal  deformation  over  time  are  also  a  source  of  star  tracker  error. 
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III.  KALMAN  FILTERING  FOR  ATTITUDE  ESTIMATION 


As  shown  in  the  preceding  section,  satellite  attitude  control  systems  estimate  the 
attitude  of  the  satellite  based  upon  measurements  from  imperfect  sensors.  As  time  passes 
the  sensors  degrade  and  the  data  they  send  to  the  attitude  control  system  degrades  as  well. 
The  Kalman  filter  provides  a  solution  to  this  problem.  In  this  section  the  Kalman  filter  is 
introduced  in  its  basic  form.  Attitude  error  representations  are  introduced  and  then  a  state 
vector  is  defined  for  attitude  estimation  purposes.  A  discrete  Kalman  filter  for  attitude 
estimation  based  upon  the  Gibbs  parameter  is  then  developed  and  the  results  of  its 
implementation  are  presented. 

A.  THE  KALMAN  FILTER 

As  alluded  to  previously,  a  Kalman  filter  is  a  recursive  mathematical  algorithm 
that  allows  one  to  estimate  the  state  of  the  system  based  upon  its  previous  states,  its 
known  dynamics,  and  knowledge  of  the  accuracies  of  the  sensors  involved.  Once 
initialized,  there  are  three  steps  to  the  filtering  process:  prediction,  measurement,  and 
correction.  A  diagram  of  these  steps  is  shown  below.  There  are  two  main  types  of 
Kalman  filters  -  regular  and  extended.  A  brief  overview  of  each  is  presented  here. 


Measurement 


* 

Initialization 

Figure  3.  Kalman  Filtering  Process  Diagram 


Correction 


17 


1.  The  Kalman  Filter 

The  Kalman  filter  may  be  used  to  predict  the  future  state  of  a  process  governed  by 
the  following  equation 

**+i  =  Fkxk  +  Guk  +  wk  (3-1) 

where  k  represents  the  time  step  and  x  is  a  vector  of  state  variables  xeH V  .  Fk  is 
an  n  x  n  matrix  known  as  the  state  transition  matrix  which  relates  the  current  state  to  the 
next  one.  The  matrix  G  is  a  nxl  matrix  that  relates  the  control  input  at  time  step  k  -  uk  - 

to  the  state  x.  The  final  term  wk  is  the  process  or  plant  noise.  This  noise  accounts  for 
random  inputs  and  inaccuracies  of  the  state  space  model  that  cause  the  plant  to  perform  in 
a  non-deterministic  manner.  The  noise  is  assumed  to  be  white  noise  with  a  nonnal 
probability  distribution 

p(w)UN(Q,Q)  (3-2) 

where  Qk  is  a  matrix  representing  the  plant  noise  w  defined  by 

Qk=cox{wk)  (3-3) 

Although  the  state  vector  represents  all  the  variables  being  estimated,  not  all  the 
estimated  variables  are  required  to  be  measured.  The  measurement  is  related  to  the  state 
by  the  following  equation 

zk=Hkxk+vk  (3-4) 

where  Hk  is  an  m  x  n  matrix  often  referred  to  as  the  measurement  matrix  or  the 

measurement  sensitivity  matrix  in  some  cases.  The  second  term  vk  represents  the 

measurement  noise  -  or  sensor  inaccuracies.  As  with  the  plant  noise,  it  is  assumed  to  be 
white  noise  with  a  normal  distribution  given  by 

p(vk)-  N(0,Rk)  (3-5) 
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where  Rk  is  the  measurement  error  covariance  matrix  defined  by 

Rk  =  covOJ  (3-6) 

It  is  important  to  note  that  though  they  share  a  similar  probability  distribution,  the 
plant  and  measurement  noise  are  completely  independent  of  one  another.  If  a  system 
under  consideration  can  be  described  by  equations  3-1  and  3-3,  then  a  discrete  Kalman 
filtering  algorithm  may  be  used  to  estimate  its  future  states. 

a.  Prediction 

The  first  phase  of  a  filtering  algorithm  is  the  prediction  phase.  In  this 
phase  at  time  step  k  a  prediction  xk+l  is  made  for  the  next  time  step  based  upon  the 
current  state  using  the  system  model.  This  prediction  is  given  by 

x~k+1=Fkxk+Guk  (3-7) 

the  minus  in  the  superscript  of  x  denotes  that  the  estimate  has  been  made 
prior  to  taking  a  measurement. 

The  error  covariance  is  also  predicted  during  the  prediction  stage.  The 
covariance  matrix  P  is  a  n  x  n  matrix  representing  the  estimate  error  covariance.  The 
covariance  matrix  is  predicted  via 

Pk'+l=FkPkFkT+Qk  0-8) 

b.  Measurement 

Once  a  prediction  for  the  next  state  has  been  made,  an  actual  measurement 
zk  is  taken.  Note  once  again  that  it  is  not  necessary  to  measure  every  state  being 
estimated.  A  quantity  known  as  the  residual  zk  is  then  calculated  by 

zk  =  zk  -Hkxk  (3-9) 

A  quantity  called  the  Kalman  gain  is  then  defined  according  to 
Kk  =PkHTk(HkPkHTk  +Rkyl  (3-10) 
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where 


pk  =E\(xk-xk)(xk-xk)T] 

and 

(3-H) 

xk  =%lv.] 

xk=E\-xk  \zk] 

c.  Correction 

The  Kalman  gain  in  conjunction  with  the  measurement  residual  are  then 
used  to  find  the  estimate  xk  via 


xk=xk+  Krsk  (3-12) 

The  Kalman  gain  is  also  used  to  update  the  covariance  matrix  by 
Pk:  =(I-KkHk)Pk7(I-KkHk))T  +  KkRkKTk  (3-13) 

Once  this  has  been  completed,  the  process  starts  over  again  -  hence  the 

tenn  recursive. 

d.  Filter  Initialization 

In  order  to  begin  the  Kalman  filtering  process  an  initial  estimate  of  both 
the  state  and  the  covariance  matrix  are  required  to  begin  the  recursion.  The  filter  can  be 
initialized  in  several  different  ways,  but  the  most  common  way  to  do  so  is  via  a  simple 
guess.  The  closer  to  the  actual  state  the  initial  guess  is,  the  faster  the  estimate  generated 
by  the  filter  will  reach  the  actual  value.  On  occasion  the  filter  is  initialized  using  the  first 
few  measurements,  but  since  it  is  easier  to  simply  make  an  educated  guess  this  method  is 
not  often  used. 

2.  The  Extended  Kalman  Filter 

Not  all  processes  that  are  desired  to  be  estimated  can  be  modeled  in  the  linear 
form  required  for  the  Kalman  filter.  For  these  processes,  the  Extended  Kalman  Filter 
(EKF)  can  be  used.  The  EKF  is  simply  a  regular  Kalman  filter  that  linearizes  about  the 
estimate  and  covariance. 

For  the  EKF,  the  state  vector  x  is  governed  by  a  non-linear  function / dependent 
upon  the  current  state,  a  control  input,  and  plant  noise 
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Xk+ 1  =f(xk,uk,wk) 


(3-14) 


with  the  measurement  z  being  expressed  as  a  non-linear  function  of  the  state,  as 

zk=h(xk,vk)  (3-15) 

Both  noise  terms  wk  and  vk  have  the  same  properties  as  with  the  linear  Kalman 

filter.  Since  it  is  impossible  to  know  the  noise  at  any  given  time  and  both  the  plant  and 
measurement  noise  are  characterized  as  white  noise,  both  the  state  and  measurement 
vectors  are  approximated  by  setting  the  noise  equal  to  zero 

xk+i  =  f(xk ,  Uk  ,0)  zk=  h(xk  ,0)  (3-16) 

Linearizing  about  the  zero  mean  non-linear  functions  in  equation  3-16  will  yield  a 
close  approximation  to  the  actual  value  of  the  function  itself  -  if  the  noise  is  truly  zero- 
mean  in  nature,  its  expected  value  must  be  zero  by  definition.  We  linearize  the  system  by 
taking  the  Jacobians  of f  and  h  with  respect  to  the  state  vector  xk ,  and  the  following 
matrices  are  obtained 


2*3 

ii 

(3-17) 

TT  _8hi 

iJ  dx,. 

(3-18) 

With  the  above  matrices  defined  the  EKF  algorithm  proceeds  in  exactly  the  same 
manner  as  the  regular  Kalman  filter. 

a.  Prediction 

The  prediction  for  the  EKF  is  done  via  the  non-linear  function/ described 
previously.  For  the  predicted  state  and  covariance  the  equations  are 

x;+l=f(xk,uk,  0)  (3-19) 

P^=FkPkFkT+Qk  (3-20) 
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b.  Measurement 

Once  a  prediction  has  been  made  a  measurement  zk  is  taken  and  as 
before,  a  residual  zk  is  calculated,  this  time  via  the  non-linear  function  h 

zk  =zk  -h(xk,0)  (3-21) 

The  Kalman  gain  is  calculated  as  before  in  equation  3-9,  with  the 
exception  that  when  using  the  EKF  H  is  the  Jacobian  of  the  non-linear  function  h  vice  the 
sensitivity  matrix  H  used  in  the  regular  Kalman  filter. 

c.  Correction  and  Initialization 

As  with  the  regular  Kalman  filter,  the  EKF  corrects  the  estimate  and  the 
covariance  matrix  based  upon  the  measurement  and  the  Kalman  gain.  The  equations 
used  to  do  this  are  identical  to  the  regular  Kalman  filter.  Also  like  the  regular  Kalman 
filter,  an  EKF  must  be  initialized  with  a  starting  value  -  like  the  correction  step  of  the 
filtering  process,  the  EKF  method  for  initialization  is  identical  to  that  used  for  the  regular 
Kalman  filter. 

B.  ATTITUDE  ERROR  REPRESENTATIONS 

When  using  Kalman  filters  for  spacecraft  attitude  estimation,  the  state  vectors 
most  frequently  consist  of  bias  and  attitude  errors.  Prior  to  developing  a  Kalman  filter  for 
attitude  estimation,  it  is  appropriate  to  cover  how  the  attitude  error  may  be  represented. 
The  two  attitude  error  representations  relevant  to  the  discussion  here  are  the  quaternion 
error  and  the  Gibbs  parameter  error  representations.  Emphasis  will  be  placed  upon  the 
latter  of  the  two. 

1.  Quaternion  Error 

If  two  reference  frames  are  slightly  offset  from  one  another,  the  error  quaternion 
Sq  represents  the  rotation  that  will  align  one  frame  with  the  other.  For  example,  if  the 
estimated  attitude  of  a  spacecraft  is  given  by  qref  and  the  measured  attitude  of  the 
spacecraft  by  a  perfect  sensor  was  q  then  Sq  would  represent  the  rotation  from  qref ,  the 
estimated  attitude,  to  q  ,  the  actual  attitude.  For  all  rotations  it  can  be  shown  that 

q  =  Sq®qref  (3-22) 
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The  ease  of  processing  the  above  equation  makes  the  error  quaternion  a  useful 
representation  for  spacecraft  attitude  error.  However,  from  a  Kalman  filtering  standpoint, 
the  error  quaternion  presents  some  problems.  Recall  that  when  defining  a  residual,  zk 

the  estimated  value  is  subtracted  from  the  measured  value.  For  illustrative  purposes, 
define  p  as  the  measured  error  quaternion  and  q  as  the  predicted  error  quaternion  based 
upon  the  current  state.  Subtracting  to  obtain  the  residual  yields 

Pi  1  kl  \Pi~q  i 

Pi  q2  p2-q2 

z=p~q= 

p2  q2  p3-q3 

_Pa\  L?4J  lP<~q< 

While  this  seems  straightforward  enough,  recall  that  attitude  quaternions  reside  on 
the  surface  of  a  hyper-sphere  defined  previously  and  shown  again  here 

q{  +  q~2  +  q\  +  ql  =  1  (2-5) 

Though  p  and  q  both  satisfy  the  requirements  of  equation  2-5  their  difference  will 
not  do  so,  and  therefore  is  not  an  attitude  quaternion.  This  becomes  a  problem  when 
using  quaternions  in  a  Kalman  filtering  algorithm,  so  converting  the  quaternion  to 
different  representations  becomes  necessary  -  as  will  be  shown  later  in  this  chapter. 

2.  Gibbs  Error 

Although  the  Gibbs  parameter  is  closely  related  to  the  quaternion  -  as  one  can 
easily  see  from  its  definition  in  chapter  two  -  it  does  not  have  the  normalization 
constraint.  The  Gibbs  parameter  for  a  particular  quaternion  resides  in  a  plane  tangent  to 
the  surface  of  the  hyper-sphere  on  which  quaternions  reside.  A  direct  mapping  between 
the  Gibbs  parameter  plane  and  the  quaternion  hyper-sphere  exists  and  is  given  by 


i 

Sq(ag)  =  (4  +  ag)  2 


(3-24) 


It  can  further  be  shown  that  the  direction  cosine  matrix  associated  with  a 
particular  5q(a)  may  be  approximated  by 

A(Sq(a ))  «  /3x3  -[ax]-^(a2IM  - aa  )  (3-25) 
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This  direct  relationship  between  the  Gibbs  parameter  and  the  quaternion  will  be 
used  in  the  next  section  to  formulate  a  Kalman  filter  for  spacecraft  attitude  estimation 
[Markley]. 

C.  FORMULATION  OF  A  GIBBS  PARAMETER  BASED  KALMAN  FILTER 

With  the  Kalman  filter  and  attitude  error  representation  introduced,  the 
formulation  of  the  Kalman  filter  may  be  developed.  The  Kalman  filter  developed  in  this 
section  estimates  the  spacecraft  attitude  error  and  the  gyroscopic  bias.  It  uses  a  six 
element  state  vector  of  the  form 


a(t ) 
bit) 


(3-26) 


where  a(t)  represents  the  attitude  error  in  Gibbs  parameters  and  b(t)  is  the 
gyroscopic  bias  error. 

1,  Filter  Overview 

The  Kalman  filter  developed  here  is  based  on  one  developed  by  F.  Landis 
Markley  [Markley].  It  treats  the  actual  attitude  as  the  quaternion  product  of  the  estimated 
quaternion  error  and  the  previous  estimate,  as  shown  here 

q(t)  =  Sq(a(t))  0  qref  (3-27) 


where  qref  is  the  attitude  estimate  from  the  previous  time  step  and  Sq(a(t)) 
represents  the  rotation  from  the  last  time  step  to  the  attitude  at  the  current  time  step. 

The  steps  for  the  filter  are  the  same  as  for  any  other  -  prediction,  measurement 
update,  and  correction.  What  is  unique  about  this  filter  is  that  the  propagation  of  the  state 
a(t) 


vector  x  = 


bit ) 


is  relatively  trivial  in  nature,  as  Markley  himself  states.  Prior  to 


developing  the  filter  itself,  it  is  of  interest  to  show  why  this  is  the  case. 

In  the  prediction  step  of  this  filter,  the  attitude  quaternion  from  the  previous  time 
step  q~ef  is  considered  the  optimal  estimate  for  the  attitude  of  the  spacecraft,  i.e.  the 

attitude  error  estimate  a~(t )  at  the  beginning  of  every  time  step  is  zero  by  definition. 
Once  a  measurement  is  taken,  q~ef  is  no  longer  the  optimal  estimate  for  the  spacecraft 
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attitude.  During  the  correction  step,  the  corrected  error  estimate  a+(t )  is  used  to 
propagate  the  old  estimate,  qri;f  so  that  it  becomes  the  new  optimal  estimate  for  the 
spacecraft  attitude  via 

q+ref  =  Sq{a+{t))  0  qKf  (3-28) 

Once  the  new  optimal  estimate  has  been  generated,  a(t )  is  immediately  reset  to 
zero.  When  the  next  step  in  the  recursion  begins,  the  new  qref  then  becomes  the  optimal 
estimate 

q~ref  =  8q(a(t))  0  qKf  =  8q{  0)  0  qKf  =  q~ref  (3-29) 

Because  the  filter  is  designed  in  this  manner,  the  actual  propagation  of  the  state 
vector  is  trivial  because  its  first  three  elements  are  always  reset  to  zero.  The  value  of 
q~ef  immediately  after  the  correction  step  of  the  filtering  process  is  the  actual  attitude 

estimate  -  the  state  vector  merely  tracks  the  error. 

2.  System  Dynamics  -  Developing  the  State  Transition  Matrix 

The  first  step  in  formulating  any  Kalman  filter  is  mathematically  modeling  how 
the  system  should  behave  and  subsequently  constructing  a  state  transition  matrix  based 
upon  this  mathematical  model.  For  a  satellite,  we  have  its  current  state  quaternion  and  its 
angular  rate  co  .  From  equation  2-12,  it  is  known  that  the  time  derivative  of  the  quaternion 
is  based  upon  the  angular  rate  of  the  body  it  is  describing. 

Placing  this  in  tenns  of  the  filter  formulation 

U  =  \  ®4V  (3-30) 

where  cbref  is  the  angular  rate  of  the  satellite  at  the  reference  attitude,  qref . 

The  filter  itself  is  not  estimating  the  quaternion  components  of  the  attitude  -  it 
estimates  the  error  in  tenns  of  Gibbs  parameters.  The  equation  for  the  propagation  of  the 
Gibbs  parameter  error  is  given  by 

1  T  1 

as  =  (73x3  +-agag\o)-(oref)--{(o+  (oref  )xag=f(a,t ) 
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(3-31) 


as  shown  by  Markley  [Markley]. 

We  now  propagate  the  estimate  from  time  nT  to  time  nT  +  T .  Therefore 

qref  =  q{nT)  =  constant  (3-32) 

and  coref  =  0  since  qref  is  constant.  Therefore,  from  equation  3-27  we  obtain 
q(nT  +  T)  =  Sq(a(nT  +  T))  0  q(nT)  (3-33) 

Propagating  in  a  discrete  manner,  ag(nT  +  T )  may  be  written  as 
ag(nT  +  T)  =  ag(nT)  + agT  (3-34) 

When  this  approach  is  taken,  equation  3-31  may  be  rewritten  as 

1  T  1 

ag=(I^+~agag)co--coxag  (3-35) 


Linearizing  about  ag  yields  the  following  equation 


a 


g 


1 

CO  —  &>x<7 
2  g 


(3-36) 


Recalling  from  equation  2-16  that  the  angular  velocity  is  actually  a  function  of  the 
measured  velocity  and  the  gyroscopic  bias  and  substituting  that  into  equation  3-36 
(assuming  a  zero  noise  component  for  the  time  being)  gives 


a  =  co  - b  —  —  (co  -b)xa 

g  meas  ^  ^  m  '  S 

—  co  —  b - — Q  a 

meas  ^  m  § 


(3-37) 


Introducing  the  term  Qm  which  is  defined  as  the  skew  of  the  measured  angular 
velocity 


Q  =  \oo  x  = 

m  L  meas 


-co.,. 


Mmeas,  0 


co,„ 


—co,„ 


—co,„ 


a 


(3-38) 


We  propagate  the  attitude  from  time  nT  to  time  nT+T  equation  3-34  may  be 
rewritten  as 


ag(nT  +  T)  =  a(nT)  +  (o)meas(nT)-b(nT)-^  Qm(nT)a(nT))T 
=  (I3a~nm(nT)T)a(nT)  +  o>mm(nT)T-b(nT)T 


(3-39) 


Assuming  a  constant  bias  over  short  periods  of  time  T,  the  state  equation  3-39 
may  be  used  to  formulate  the  state  transition  matrix  for  the  Kalman  filter 


4  = 


/  -  —  Q  T  -I  T 

1 3x3  2  3*3 


0 


3x3 


I 


3x3 


(3-40) 


which  will  propagate  the  state  as 

x(nT  +  T)  =  Adx(nT)  +  u(nT)  +  w(nT)  (3-41) 

This  is  in  the  standard  state  space  form  so  we  are  able  to  apply  standard  Kalman 
filtering  algorithms  with  the  control  input  modeled  as 

u  =  ooT  (3-42) 

3.  System  Sensors  -  Developing  the  Measurement  Sensitivity  Matrix 

With  the  state  transition  matrix  developed  from  the  system  dynamics,  it  becomes 
necessary  to  determine  how  the  filter  will  take  an  external  measurement  and  incorporate 
it  into  the  filtering  algorithm  via  the  measurement  sensitivity  matrix.  The  term  which 
drives  the  ‘correction’  equation  is  the  error  between  the  actual  measurement  and  the 
predicted  measurement.  In  particular  in  our  problem  we  have  the  actual  orientation  'v 
from  the  star  trackers  library  in  the  inertial  coordinate  frame.  At  the  same  time,  given  an 
estimate  of  the  spacecraft  orientation  q  and  the  observed  star  hv  from  the  star  tracker  we 
can  compute  and  estimate 

'v  =  A(q)bv  (3-43) 

In  particular  let 

q  =  Sq(a)®  q+ (nT)  (3-44) 
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so  the  error  'v  -  'v  can  be  written  as 

'v-'v=  'v-  A(Sq(a))A(q+(nT))bv  (3-45) 

Substituting  equation  3-35  into  the  above  equation  produces  a  relationship 
between  the  Gibbs  error  a,  qref ,  hv  and  'v 

*v  =  [4,3  -[ax]-I(a2/3x3 -aar)]A(qrefyv  (3-46) 

In  the  above  equation  qref  is  the  estimate  from  the  previous  time  step. 

Conducting  a  first  order  expansion  with  respect  to  the  Gibbs  error  about  bv  yields  the 
following  equation 

p,h  _ 

Kbv )  =  h(bv )  +  —  [ax]' V  =  h{bv )  +  Haa  (3-47) 

d  v 


From  this,  Ha  is  shown  to  be 

Ha  =  [*Fx]  (3-48) 

d  v  »= 

With  the  above  relations  developed,  the  measurement  sensitivity  matrix  can  be 
formed  by  relating  them  to  the  state  vector  x 


H  =  \Ha0  3x3]  (3-49) 

4.  The  Filter  Equations 

With  the  state  transition  and  measurement  sensitivity  equations  constructed,  it  is 
now  possible  to  summarize  and  implement  the  Kalman  filter.  In  this  version  of  the  filter, 
the  prediction  and  correction  steps  are  combined  for  expediency. 

The  Kalman  filter  is  built  on  the  following  state  space  model,  developed  above 


x(n  + 1)  =  2  m  3x3  x(n)  +  +  w  =  Ax  +  Bu  +  w 

0  0  L01x3. 

L  *-^3x3  ^3x3  J 

"  fit, 

y(n)=  -TbZ  [/,vx]03x3]  x(n)  +  v  =  Hx  +  v 

O  V  b= 


(3-50) 
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where  x(n )  = 


and  u(n)  =  (omeJ  . 


a 

b 


The  following  equations  are  used  to  estimate  the  state  and  the  covariance 


x. 


k+ 1 


- Q  T  -U.J 


3x3  ~  m 


0 


3x3 


3x3 


‘3x3 


xk  +K(zk  -[[  ’vx]  03 x3]xk)  +  u 


(3-51) 


=  Adxk  +K(zk-Hxk)  +  u 


Pk+l  =  AAK  +  Q-Kk{HPkHT+R)KkT  (3-52) 


Note  that  the  above  two  equations  differ  slightly  from  the  Kalman  filtering 
equations  derived  earlier.  For  implementation  purposes,  the  prediction  and  correction 
steps  were  combined.  It  is  also  important  to  note  here  that  the  measurement  zk  is  the 
measured  error,  which  is  given  by 


z 


k 


=  h V-- 


(3-53) 


The  Kalman  gain  K  used  in  equation  3-60  is  calculated  by 

Kk=AdPkHT(HPkHT+RTx  (3-54) 

A  diagram  of  the  filtering  process  is  shown  in  Figure  4.  The  reference  quaternion 
which  is  updated  after  the  filter  estimates  the  attitude  error  gives  the  optimal  estimate  of 
the  attitude  of  the  spacecraft  after  it  is  updated  with  the  information  from  the  filter.  Once 
the  reference  quaternion  becomes  the  optimal  estimate,  the  attitude  error  estimate  is  reset 
to  zero. 
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Figure  4.  Gibbs  Parameter  Kalman  Filter  Flow  Diagram 


5.  Testing  the  Filter 

Once  developed,  the  filter  was  implemented  in  MATLAB  and  tested  using  both 
‘clean’  and  corrupted  quaternion  measurements  taken  from  a  simulated  satellite  that 
drives  itself  to  nadir  pointing  after  an  initial  offset.  The  gyroscopes  in  this  simulation  had 
no  bias.  Clean  quaternions  and  gyroscope  readings  were  used  on  this  initial  run  to  see  if 
the  Kalman  filter  would  function  correctly.  The  results  are  shown  in  figures  five  and  six. 
Though  some  quaternion  error  is  present,  it  is  on  the  order  of  1 0  6  and  is  attributed  to 
numerical  error  only.  The  estimated  bias  is  of  the  same  order  and  is  also  attributed  to 
numerical  errors  and  the  fact  that  a  non-zero  value  for  the  bias  was  used  to  initialize  the 
filter.  Because  of  the  low  level  of  error,  it  was  concluded  that  the  filter  functions 
correctly. 
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Quaternion  Errors 
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Figure  5.  Quaternion  Errors  For  Gibbs  Parameter  Kalman  Filter  With  Zero  Bias  and  No 

Measurement  Noise 
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Figure  6.  Bias  Estimated  By  Gibbs  Parameter  Kalman  Filter  with  Zero  Bias  and  No 

Measurement  Noise 
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IV.  STAR  GAP  ERROR  MITIGATION 


As  a  satellite  travels  in  its  orbit  about  the  Earth,  there  will  be  times  when  none  of 
its  star  trackers  are  able  to  feed  an  attitude  quaternion  into  the  attitude  control  system. 
There  are  several  possible  reasons  for  these  star  gaps,  including  solar  and  lunar  entry  into 
the  field  of  view  (FOV)  of  the  star  tracker,  processing  time  for  the  tracker  to  find  and 
lock  on  a  new  star(s)  once  the  previous  one(s)  exit  the  FOV,  or  due  to  a  reorientation 
maneuver  by  the  satellite  -  among  others.  These  gaps  have  marked  effects  upon  attitude 
estimation  algorithms.  In  this  section  the  effects  of  star  gaps  are  examined  -  both  for  the 
Gibbs  parameter  based  Kalman  filter  developed  in  Chapter  Three  and  an  Euler  angle 
based  Kalman  filter  developed  by  Palermo  using  the  same  data  [Palermo].  Methods  to 
mitigate  the  effects  of  these  star  gaps  are  presented  along  with  simulated  results  of  their 
implementation. 

A.  EFFECTS  OF  STAR  GAPS 

When  a  star  gap  occurs,  the  effect  it  has  on  the  attitude  estimation  algorithm  is 
highly  dependent  upon  the  formation  of  that  algorithm.  Regardless  of  the  type  of 
algorithm  in  use,  a  longer  star  gap  equates  to  a  larger  error  as  will  be  shown  in  the  results 
later  in  this  section. 

1.  Euler  Angle  Based  Kalman  Filter 

Prior  to  examining  how  star  gaps  affect  the  Euler  angle  based  Kalman  filter  it  is 
useful  to  briefly  examine  its  formulation.  For  a  full  derivation  of  the  equations,  the 
reader  is  referred  to  [Palermo].  When  a  star  gap  occurs  in  the  Euler  angle  based  Kalman 
filter,  the  attitude  estimate  diverges  for  all  four  elements  of  the  quaternion.  Upon 
reacquisition  of  a  star,  there  is  a  brief  spike  in  the  attitude  error  then  the  measurement 
updates  drive  the  error  back  down  -  eventually.  The  results  of  a  200  second  star  gap  on 
an  attitude  estimator  based  upon  Euler  angles  are  shown  in  the  following  figure.  The 
estimator  is  tracking  noisy  measurements  for  both  attitude  and  angular  rate. 
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Quaternion  Error 


Figure  7.  Euler  Angle  Based  Kalman  Filter  with  200  Second  Star  Gap 

2.  Gibbs  Parameter  Based  Kalman  Filter 

The  Gibbs  Parameter  based  Kalman  Filter  behaves  in  a  very  different  manner. 
When  a  star  gap  occurs,  the  estimator  is  able  to  continue  the  track  with  excellent  accuracy 
because  it  is  still  receiving  angular  rate  data.  The  measured  angular  rate  data  allows  the 
estimator  to  propagate  the  reference  quaternion  qref  with  little  loss  of  accuracy  due  to  the 

star  gap.  The  graph  was  made  using  the  exact  same  data  for  a  200  second  star  gap  as  was 
used  for  the  Euler  angle  based  Kalman  filter.  Figure  8  shows  the  quaternion  errors  for 
the  same  track. 
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Figure  8.  Estimated  Quaternion  Elements  vs.  Actual  for  Gibbs  Parameter  Based  Kalman 

Filter  (Estimate  Shown  in  Blue)  with  200  Second  Star  Gap 
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Figure  9.  Quaternion  Error  for  Gibbs  Parameter  Based  Kalman  Filter  with  200  Second  Star 

Gap 
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B.  ADAPTIVE  COVARIANCE 

One  technique  for  mitigating  the  effects  of  a  star  gap  is  adapting  the  different 
covariance  matrices  intrinsic  to  the  Kalman  filtering  process.  The  three  matrices  dealing 
with  covariance  are  the  filter  covariance  matrix  P,  the  plant  covariance  matrix  Q,  and  the 
measurement  covariance  matrix  R.  The  P  matrix  is  predicted  and  corrected  based  upon 
measured  values  and  the  Kalman  gain  -  adapting  it  based  upon  situational  parameters 
would  significantly  change  the  way  in  which  the  filter  propagates  the  estimate.  The  plant 
and  measurement  covariance  matrices,  however,  remain  constant  throughout  the  filtering 
process.  The  performance  of  any  Kalman  filter  may  be  modified  by  changing  their  value. 
By  doing  so  adaptively  based  upon  external  conditions,  it  is  possible  to  mitigate  some  of 
the  attitude  error  caused  by  intennittent  star  gaps.  Adaptation  of  the  plant  and 
measurement  covariance  are  examined  individually  and  then  examined  in  conjunction 
with  one  another. 

C.  THE  PLANT  COVARIANCE  MATRIX 

The  plant  covariance  Q  may  be  thought  of  as  a  measure  of  how  well  the  model 
emulates  the  actual  dynamics  of  the  quantities  being  estimated.  The  more  precise  the 
model,  the  lower  the  values  contained  in  the  Q  matrix.  A  star  gap  might  be 
conceptualized  as  a  system  modeled  in  an  extremely  poor  manner,  i.e.  a  system  with  an 
extremely  large  plant  noise.  Using  this  conceptualization,  a  trigger  was  developed  that 
increased  the  plant  covariance  when  a  star  gap  occurred.  This  covariance  trigger  was 
implemented  on  both  the  Euler  and  Gibbs  parameter  based  Kalman  Filters. 

1.  Adapting  Plant  Noise  for  the  Euler  Angle  Based  Kalman  Filter 

The  Euler  angle  based  Kalman  Filter  used  for  testing  the  adaptive  covariance 
concept  was  developed  by  Palenno  using  SIMULINK  [Palermo].  This  fonnat  was 
maintained,  but  modified  to  adapt  the  covariance  matrices.  The  following  figure  shows 
the  trigger  that  was  developed  for  this  purpose.  The  figure  shown  is  the  measurement 
covariance  trigger  -  the  plant  noise  covariance  trigger  is  identical.  An  initial  matrix 
serves  as  the  covariance  trigger  until  a  certain  point  -  user  determined  -  at  which  it  is 
increased  by  a  user  determined  factor.  After  a  certain  time  period  has  elapsed,  the 
covariance  matrix  is  reset  to  its  initial  value.  A  digital  clock  was  used  to  do  this,  as  can 
be  seen  in  the  diagram. 
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Figure  10.  Measurement  Covariance  Trigger 

Once  the  covariance  trigger  was  implemented,  it  was  tested  using  several  different 
gains  for  the  200  second  star  gap  under  consideration.  Note  that  the  trigger  was  timed  to 
increase  the  covariance  just  before  the  star  gap  and  decrease  it  immediately  just  after  the 
star  gap  -  assuming  that  the  measurements  preceding  the  star  gap  were  stored  in  memory 
and  recycled  through  with  the  new  value  of  Q  once  the  star  gap  was  detected.  Both  the 
magnitude  and  the  directions  of  the  rotations  varied  for  each  gain  used.  The  gains  for  the 
plant  noise  covariance  matrix  were  increased  logarithmically  for  the  duration  of  the  star 
gap.  The  following  figures  shows  the  mean  quaternion  errors  plotted  against  the  gains 
applied  to  the  plant  noise  matrix  Q  for  the  duration  of  the  star  gap. 
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Mean  Error  for  Quaternion  Elements  with  Plant  Noise  Injection  During  Star  Gap 
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Figure  1 1 .  Mean  Quaternion  Error  for  Plant  Noise  Injection  During  Star  Gaps 

2.  Adapting  Plant  Noise  for  the  Gibbs  Parameter  Based  Kalman  Filter 

The  Kalman  filter  based  on  Gibbs  Parameters  was  built  in  MATLAB  without 
using  SIMULINK,  making  any  adaptation  of  the  Q  matrix  easy  to  accomplish.  Unlike 
the  Euler  angle  based  Kalman  filter,  adapting  the  plant  noise  with  the  Gibbs  Parameter 
based  Kalman  filter  produced  no  noticeable  effect. 

D.  THE  MEASUREMENT  COVARIANCE  MATRIX 

The  measurement  covariance  matrix  R  provides  the  filter  with  the  known 
accuracies  of  the  sensors  being  used  to  take  the  measurements  z.  Conceptually,  it  follows 
that  when  a  sensor  is  producing  no  measurement  its  function  should  come  into  question. 
This  was  done  for  both  types  of  Kalman  filters  discussed  in  this  paper. 

1.  Perturbation  of  Measurement  Noise  Covariance  Matrix  for  an  Euler 
Angle  Based  Kalman  Filter 

Adapting  the  measurement  covariance  matrix  accordingly  using  the  trigger  shown 
in  figure  eight  produced  the  following  results  when  applied  to  the  Euler  angle  based 
Kalman  filter. 
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x  -jo  4  Mean  Error  for  Quaternion  Elements  with  Measurement  Noise  Injection  During  Star  Gap 


Figure  12.  Mean  Quaternion  Error  for  Measurement  Matrix  Adaptation  During  Star  Gap 

2.  Perturbation  of  Measurement  Noise  Covariance  Matrix  for  a  Gibbs 
Parameter  Based  Kalman  Filter 

As  with  the  plant  noise  injection,  measurement  noise  injection  for  the  Gibbs 
parameter  based  Kalman  filter  produced  no  difference  with  regards  to  estimation  error 
encountered  during  a  star  gap. 

E.  SIMULTANEOUS  ADAPTATION  OF  COVARIANCE  MATRICES 

If  adjusting  the  sensitivity  of  the  filter  to  plant  and  measurement  noise 
individually  affects  the  average  quaternion  error  for  a  Kalman  filter,  it  follows  that 
adapting  both  simultaneously  would  do  so  as  well.  This  was  accomplished  via  use  of  two 
triggers  operating  in  parallel,  as  shown  previously  in  figure  eight.  The  results  are  shown 
in  the  following  figure  -  a  significant  improvement  over  adapting  either  matrix 
individually. 
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Average  Quaternion  Error  with  Plant  and  Measurement  Noise  Injection 
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Figure  13.  Average  Quaternion  Error  with  Simultaneous  Plant  and  Measurement  Covariance 

Matrix  Adaptation 

Adaptation  of  the  plant  noise  and  measurement  noise  matrices  is  an  option  for 
mitigating  attitude  estimation  error  in  some  cases,  as  has  been  shown  in  the  preceding 
sections.  While  it  does  not  eliminate  the  error  induced  due  to  star  gaps,  it  does  have  the 
potential  to  reduce  the  amount  of  error  encountered. 
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V.  DYNAMIC  GYROSCOPE  AND  RATE  GYRO  UPSETS 


Mechanical  rate  gyroscopes  have  historically  been  a  weak  point  in  satellite 
design.  Their  high  failure  rate  coupled  with  their  noisy  output  and  importance  in  attitude 
estimation  makes  them  a  liability  from  a  reliability  standpoint.  Improved  gyroscopes 
have  been  developed  -  such  as  the  laser  gyroscope  -  but  many  systems  currently  in  orbit 
are  still  relying  on  mechanical  gyroscopes  which  are  both  noisy  and  prone  to  failure. 

This  chapter  introduces  the  dynamic  gyroscope  -  a  software  alternative  to  the  rate 
gyroscope,  shows  how  rate  gyroscope  upsets  can  affect  attitude  estimation  algorithms, 
and  then  integrates  the  dynamic  gyroscope  with  attitude  estimation  algorithms  to  mitigate 
the  effects  of  rate  gyroscope  upsets. 

A.  THE  RATE  GYROSCOPE  UPSET 

The  rate  gyroscope  upset  is  a  phenomenon  that  has  been  occurring  on  certain 
satellites  in  orbit.  At  certain  random  points  in  time  the  rate  gyroscope  will  send  a 
measurement  to  the  attitude  control  system  that  is  100  to  1000  times  the  actual  reading. 
This  wreaks  havoc  on  the  control  system,  as  can  be  seen  in  the  following  two  figures. 
Both  figures  were  produced  with  the  same  data  from  a  simulated  satellite  -  with  a 
simulated  rate  gyroscope  upset.  The  simulated  rate  gyroscope  upset  was  of  two  seconds 
duration  with  a  measurement  1000  times  the  correct  reading  being  sent  to  the  attitude 
estimation  system.  As  can  be  seen  by  looking  at  the  figures,  the  Gibbs  Parameter  based 
Kalman  Filter  developed  earlier  performs  extremely  well  as  compared  to  the  Euler  angle 
based  filter.  After  a  brief  jump  in  estimator  error,  the  attitude  error  for  the  Gibbs 
Parameter  Kalman  filter  returns  to  the  order  oflO4  while  the  Euler  Angle  Kalman  filter  is 
unable  to  recover  any  semblance  of  accuracy  after  an  upset  occurs. 
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Figure  14. 


Figure  15. 


Quaternion  Error  -  Euler  Angle  Model  with  Two  Second  Rate  Gyroscope  Upset 
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Quaternion  Error  with  2  second  Rate  Gyroscope  Upset  -  Gibbs  Parameter  Based 

Kalman  Filter 


42 


Rate  gyroscope  upsets  of  longer  and  longer  lengths  were  tested  on  the  Gibbs 
Parameter  Kalman  filter,  as  shown  in  figure  16  and  up  to  a  certain  point  the  attitude 
estimator  was  able  to  recover  its  pre -upset  estimation  accuracy.  After  a  certain  time 
period  the  estimator  does  lose  its  ability  to  recover  from  an  upset,  as  shown  in  figure  17 
when  the  rate  gyroscope  upset  was  150  seconds  in  length.  However,  its  superiority  over 
the  Euler  Angle  Kalman  filter  remains  apparent. 


Quaternion  Error  -  Gibbs  Parameter  Model  with  95  Second  Rate  Gyroscope  Upset 


Figure  16.  Gibbs  Parameter  Based  Kalman  Filter  with  95  second  Rate  Gyroscope  Upset 


43 


Quaternion  Error  -  Gibbs  Parameter  Model  with  100  Second  Rate  Gyroscope  Upset 


Figure  17.  Gibbs  Parameter  Based  Kalman  Filter  with  a  150  Second  Rate  Gyroscope  Upset 

A  bit  of  analysis  revealed  that  the  covariance  matrix  for  the  Kalman  Filter  was 
positive  semi-definite  for  all  in  which  it  was  able  to  recover  from  the  rate  gyroscope 
upset.  If  the  covariance  matrix  became  non-positive  semi-definite  during  the  rate 
gyroscope  upset,  the  filter  would  not  recover. 

B.  THE  DYNAMIC  GYROSCOPE 

It  is  possible  to  mathematically  model  and  predict  all  of  the  external  torques  on  a 
satellite  with  some  degree  of  accuracy  -  some  more  than  others.  Recall  from  basic 
physics  that  the  torque  A  of  a  satellite  during  a  time  period  may  be  expressed  as 

Sdt  (5-1) 

Equation  5-1  may  further  be  broken  down  because  the  torques  on  a  satellite  may 
be  expressed  as  the  sum  of  its  external  torques  -  which  may  be  mathematically  modeled 
-  and  its  internal  torques  from  any  momentum  exchange  devices.  Equation  5-1  may 
therefore  be  written  as 
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N„m«=N„,  +  NmD  =  {'“fidl  (5-2) 

J,k 

It  follows  that 

A,  +  Nmed  =  J'“  Sdt  =  H,m  -Hh  (5-3) 

Assuming  that  any  change  in  the  moment  of  inertia  between  tk+[  and  tk  will  be 
negligible  and  recalling  that 

Hk=Ia>k  (5-4) 

It  follows  that 

-Htk  =I(a>k+1-3k) 

~^=I~\Htki-Htk)  (5-5) 

®*+i  =®k+r\Htk+i-Htt) 

Inserting  equation  5-3  into  equation  5-5  gives  a  way  to  predict  the  angular 
momentum  for  tk+l  given  only  known  or  mathematically  modeled  data 

1  =  I  ( ^ ext  +  Nmed)  (5-6) 

Using  equation  5-6  or  some  variant  thereof,  the  angular  momentum  of  a  satellite 
may  be  obtained  analytically.  An  algorithm  that  calculates  the  angular  rate  in  this 
manner  is  referred  to  as  a  dynamic  or  a  pseudo-gyroscope.  Such  a  gyroscope  was 
developed  by  Aerospace  Corporation  and  a  variant  of  this  was  implemented  by  Palenno 
[Palermo].  The  following  figure  shows  a  comparison  of  a  simulated  real  gyroscope  and 
the  Palermo  dynamic  gyroscope.  The  dynamic  gyroscope  produces  a  reading  very  close 
to  that  produced  by  the  real  gyroscope.  While  rate  error  is  not  absent,  the  amount  of  error 
-  shown  in  figure  19  -  is  such  that  the  resulting  attitude  error  will  be  manageable. 
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Real  Gyroscope  Rates 


Figure  18.  Comparison  of  Real  and  Dynamic  Gyroscope  Performance 


Dynamic  Gyroscope  -  Gyroscope  Difference 


Figure  19.  Difference  Between  Real  and  Dynamic  Gyroscope  Readings 
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Though  convenient,  dynamic  gyroscopes  are  not  free  of  problems.  Any 
inaccuracy  in  knowledge  of  the  moment  of  inertia  of  the  spacecraft  will  result  in  angular 
rate  errors,  as  can  be  seen  from  equation  5-6.  Any  expenditure  of  fuel,  movement  of  fuel 
due  to  any  slew  maneuver  or  satellite  motion,  or  any  movement  of  the  solar  arrays  - 
among  other  things  -  will  change  the  moment  of  inertia  of  a  satellite,  thereby  causing 
some  error.  System  identification  algorithms  can  be  used  to  help  mitigate  this  problem, 
but  such  algorithms  are  not  perfect  -  error  will  always  remain. 

C.  INTEGRATING  THE  DYNAMIC  GYROSCOPE  WITH  AN  ATTITUDE 

ESTIMATOR 

By  using  a  dynamic  gyroscope  in  conjunction  with  a  real  mechanical  gyroscope  it 
is  possible  to  mitigate  the  effect  of  a  rate-gyroscope  upset.  As  seen  in  the  above  figures, 
the  dynamic  gyroscope  is  able  to  produce  data  quite  close  to  that  of  a  real  gyroscope. 
Given  the  fact  that  the  Gibbs  Parameter  based  Kalman  filter  is  extremely  robust  with 
regards  to  recovering  from  rate  gyroscope  upsets  and  the  fact  that  any  mechanism  for 
switching  from  regular  gyroscope  readings  to  measurements  from  the  dynamic  gyroscope 
will  be  imperfect,  the  Gibbs  Parameter  based  Kalman  filter  is  an  ideal  estimator  with 
which  to  integrate  the  dynamic  gyroscope.  The  following  figure  shows  the  quaternion 
error  for  a  Gibbs  parameter  based  Kalman  filter  with  a  95  second  rate  gyroscope  upset  of 
1000  times  the  actual  reading.  The  dynamic  gyroscope  in  this  scenario  took  over 
immediately  so  no  upset  readings  entered  the  estimator. 
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Figure  20.  Comparison  of  Estimator  Performance  with  and  without  Dynamic  Gyroscope 

input  during  a  95  second  Rate  Gyroscope  Upset 

As  can  be  seen,  when  a  dynamic  gyroscope  is  in  use  and  no  data  from  the  real 
gyroscope  during  its  upset  period  reaches  the  estimator,  the  amount  of  error  encountered 
is  only  slightly  increased.  As  mentioned  earlier,  this  would  be  an  ideal  scenario.  It 
would  be  impossible  to  design  a  perfect  algorithm  where  the  dynamic  gyroscope  data 
would  replace  the  real  gyroscope  measurements  immediately.  Some  upset  data  would 
always  have  a  chance  of  getting  into  the  estimator.  The  robustness  of  the  Gibbs 
parameter  based  Kalman  filter  to  short  rate  gyroscope  upsets  was  shown  previously.  The 
following  figure  shows  a  scenario  in  which  the  dynamic  gyroscope  takes  over  two 
measurements  after  the  commencement  of  the  rate  gyroscope  upset.  Notice  the  minor 
perturbations  in  the  error  -  the  estimator  rapidly  converges  back  to  an  estimate  with  error 
on  the  order  of  10  4 . 
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Quaternion  Errors  -  Rate  Gyroscope  Upset  with  Imperfect  Timing  of  DG  Takeover 


Figure  2 1 .  Quaternion  Error  with  Dynamic  Gyroscope  Takeover  Occurring  Two  Seconds 

after  Rate  Gyroscope  Upset  Commencement 

By  creating  a  switching  algorithm  with  a  high  reliability,  rate  gyroscope  upsets 
may  be  experienced  by  an  attitude  estimator  using  a  Gibbs  parameter  based  Kalman  filter 
with  little  loss  of  estimator  accuracy,  as  can  be  seen  in  the  figures  above. 
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VI.  CONCLUSIONS 


The  Gibbs  parameter  formulation  of  the  Kalman  filter  provides  an  accurate 
estimate  of  spacecraft  attitude.  It  is  far  more  robust  an  estimator  when  confronted  with 
either  a  star  gap  or  rate  gyroscope  upset  than  the  Euler  angle  estimator. 

A.  SUMMARY 

A  Kalman  filter  based  upon  the  Gibbs  parameterization  of  spacecraft  attitude  is 
developed  and  analyzed.  A  comparison  between  this  formulation  and  a  Kalman  filter 
developed  by  Palermo  is  conducted  under  different  operating  scenarios.  The  dynamic 
gyroscope  is  shown  to  be  a  viable  substitute  for  mechanical  gyroscopes  during  periods  of 
rate  gyroscope  upsets.  The  Gibbs  parameter  based  Kalman  filter  used  in  conjunction 
with  the  dynamic  gyroscope  is  shown  to  be  an  excellent  method  to  mitigate  the  effects  of 
rate  gyroscope  upsets  and  produces  far  superior  results  than  the  Euler  angle  fonnulation 
of  the  Kalman  filter  when  faced  with  a  rate  gyroscope  upset. 

B.  RECOMMENDATIONS 

It  is  strongly  recommended  that  future  work  be  done  on  attitude  estimation  in  two 
specific  areas  at  a  minimum:  Kalman  Filtering  and  Unscented  Filtering. 

1.  Kalman  Filtering 

The  Kalman  filter  developed  in  this  work  was  limited  to  six  states  for  purposes  of 
the  research.  It  is  possible  -  and  somewhat  desirable  in  some  cases  -  to  have  a  much 
larger  state  vector  estimating  other  factors  impacting  spacecraft  attitude  estimation  such 
as  scale  factors  of  gyroscopes  and  star  trackers,  alignment  errors  in  the  sensor  to  body 
transfonnation  matrices,  misalignment  of  gyroscopes,  and  errors  in  the  ECI  to  body 
coordinate  transformations  -  among  others.  Grey,  Kolve,  Herman,  and  Westerlund  have 
developed  such  a  filter  that  is  used  in  on-orbit  calibrations  [Gray].  The  Aerospace 
Corporation  has  developed  a  similar  filter.  It  is  recommended  that  such  filters  and  the 
modeling  to  support  their  implementation  be  developed  and  studied  -  implementing  the 
Gibbs  parameter  formulation  for  attitude  error  estimation  into  these  models. 

2.  Unscented  Filtering 

As  mentioned  in  the  description  of  the  EKF,  it  is  based  upon  a  linearization  of 
non-linear  functions.  This  linearization  is  the  fundamental  flaw  of  the  EKF.  Julier  et  al 
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have  developed  an  alternative  to  the  EKF  based  on  the  premise  that  it  should  be  easier  to 
approximate  a  Gaussian  distribution  than  to  approximate  an  arbitrary  non-linear  function. 
Markley  and  Crassidis  have  adapted  this  filtering  method  to  spacecraft  attitude  estimation 
and  it  is  recommended  that  research  into  this  method  as  an  alternative  to  Kalman  filtering 
be  conducted  [Crassidis], 
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